2.1 Read CSV

The Semptember Arctic sea ice extent time series is accessed in CSV format directly via an FTP archive using the read.csv() function.

data <- read.csv("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/north/monthly/data/N_09_extent_v3.0.csv")
head(data)
##   year mo     data.type  region extent area
## 1 1979  9       Goddard       N   7.05 4.58
## 2 1980  9       Goddard       N   7.67 4.87
## 3 1981  9       Goddard       N   7.14 4.44
## 4 1982  9       Goddard       N   7.30 4.43
## 5 1983  9       Goddard       N   7.39 4.70
## 6 1984  9       Goddard       N   6.81 4.11

2.2 Convert Units

The extent values are converted from millions of square kilometers (x10⁶ km²) to millions of square miles (x10⁶ mi²). The result is rounded to 2 digits.

data$extent_sqmi <- round(data$extent*0.386, digits = 2)

2.3 Create Line Plot

The ggplot2 package allows fine-tuned control of all plot elements, such as colors, axis breaks, and axis titles.

pLine <- ggplot(data = data, aes(x = year, y = extent_sqmi)) +
  geom_line(color = "#1E78B4", size = 0.75) +
  geom_point(color = "#1E78B4", size = 2.5) +
  scale_y_continuous(limits = c(0,3.5),
                     expand = c(0,0),
                     breaks = seq(0, 3.5, 0.5)) +
  scale_x_continuous(limits = c(1975, 2022),
                     expand = c(0,0),
                     breaks = seq(1975, 2020, 5)) +
  labs(y = "Extent\n[million square miles]",
       title = "Average September Artic Sea Ice Exent") +
  theme(axis.title.x = element_blank(),
        panel.grid = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(),
        plot.title = element_text(hjust = 0.5, size = 16))

pLine

3.1 Retrieve Data

First, the two ancillary Natural Earth layers are imported using the ne_download() function from the rnaturalearth package.

land <- ne_download(scale = "medium", type = "land", category = "physical", returnclass = "sf")
coast <- ne_download(scale = "medium", type = "coastline", category = "physical", returnclass = "sf")

Polygon geometries of September 1979 and 2018 sea ice extent are downloaded from the Sea Ice Index FTP archive and read using the st_read function.

ice79 <- st_read("./data/extent_N_197909_polygon_v3.0.shp")
ice18 <- st_read("./data/extent_N_201809_polygon_v3.0.shp")

3.2 Transform Data

Several steps are needed to prepare the ice, land, and coast layers for mapping:

1. Extract the native projection from the ice layer and modify the longitude of projection center
2. Crop the layers to include only geometries north on 20 degrees north latitude
3. Reproject all layers to the new custom projection
4. Create a layer to represent the visible extent of the forthcoming maps (everything north of 50 degrees north latitude)
5. Clip the ice, land, and coast layers to the visible area of interest

#[1] Extract the proj4string from the ice79 layer and change the longitude of projection center to -110
nsidc_proj <- st_crs(ice79)[[2]]
nsidc_proj
## [1] "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"

nsidc_proj <- gsub("lon_0=-45", "lon_0=-110", nsidc_proj)

#[2-3] Crop the land and coasts at 20N latitiude; reproject ice, land, and coasts to nsidc_proj
land %>% st_crop(xmin = -180, xmax= 180, ymin = 20, ymax = 90) %>% st_transform(crs = nsidc_proj) -> land2 
coast %>% st_crop(xmin = -180, xmax= 180, ymin = 20, ymax = 90) %>% st_transform(crs = nsidc_proj) -> coast2
ice79 <- st_transform(ice79, crs = nsidc_proj)
ice18 <- st_transform(ice18, crs = nsidc_proj)

#[4] Draw a line from the North Pole to 50N and get the nsidc_proj length of that line
line_coords <- matrix(c(0,90,0,50), ncol=2, byrow=TRUE)
st_linestring(line_coords) %>% st_sfc(crs = 4326) %>% st_transform(crs = nsidc_proj)-> line
radius <- as.numeric(st_length(line))

#[4 continued] Create a cicle using the radius that captures everything north of 50N
st_point(x = c(0,0)) %>% st_buffer(dist = radius) %>% st_sfc(crs = nsidc_proj) -> circle

#[5] Clip the ice, land, and coast layers to the circle
land3 <- st_intersection(land2, circle)
coast3 <- st_intersection(coast2, circle)
ice79 <- st_intersection(ice79, circle)
ice18 <- st_intersection(ice18, circle)

3.3 Create Maps

A function named plot_sea_ice is witten to draw the circle, land, coast, ice, and cicle border in that order. The two years are then arranged side by side for comparison.

plot_sea_ice <- function(x, year, value, label_size, line_width){
  p <- ggplot() + geom_sf(data = circle, fill = "#A6CEE3", color = NA) +
    geom_sf(data = land3, fill = "#B2DF8A", color = NA) +
    geom_sf(data = coast3, size = line_width) +
    geom_sf(data = x, fill = "white", size = line_width) +
    geom_sf(data = circle, fill = NA, size = line_width) +
    coord_sf(datum = NA, clip = "off") +
    theme(panel.background = element_blank(),
          axis.title = element_blank(),
          plot.margin = margin(r=2.5, unit = "cm")) +
    scale_x_continuous(expand = c(0,0)) +
    scale_y_continuous(expand = c(0,0)) +
    annotate("text", x = radius, y = -radius,
             label = paste0(year, "\n", value, " million\nsquare miles"),
             hjust = 0, vjust = 0, size = label_size)
  
  return(p)
}

p79 <- plot_sea_ice(ice79, "1979", data$extent_sqmi[1], 4, line_width = 0.25)
p18 <- plot_sea_ice(ice18, "2018", tail(data$extent_sqmi, n = 1), 4, line_width = 0.25)

plot_grid(p79, p18)